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The problem of heat transfer due to laminar flow of a viscous fluid in a channel is 
studied under the assumption that there is a parabolic distribution of velocity. The effect 
of axial temperature changes are considered and the solution is based on the simpler situa- 
tion where axial effects are discussed. The solution, obtained by the method of least squares, 
is represented in terms of a set of nonorthogonal characteristic functions. These functions 
and the corresponding characteristic values are determined by numerical integration employ- 
ing the Runge-Kutta procedure. Finally, asymptotic developments are obtained which 
are useful in the limiting cases. 



1 . Statement of Problem 

The present paper is concerned with the steady 
state problem of heat transfer in a tube. The fluid 
moves with a prescribed velocity profile, which is 
parabolic in the present case. Thus the end effect is 
not considered in the velocity profile. The heat- 
transfer is small and is supposed not to influence 
the fluid motion. The wall of the tube is assumed 
to be kept at a fixed temperature, 6 U while the 
temperature of the fluid entering the tube is fixed 
to be O . The velocity distribution in the tube 
whose radius is unity is defined by 



v x =2v m (l-r 2 ), 0<r<l, 



(1.1) 



where r is the radial distance from the center and 
v x is taken along the #-axis which is the axis of 
symmetry of the tube and v m is the average value 
of v x . The equation of heat transfer under these 
conditions is 

* bx =K \ d^+r br + Wf (L2) 

with the boundary conditions 



0=0 O for £=0 
finite for x= °° 
0=0! for r=l. 



(1.3) 



This problem has been considered [l] 4 under the 
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assumption that temperature changes in the in- 
direction will be negligible so that the term d 2 6/dx 2 
may be omitted. However, recent studies [2,3], 
which have included the effects of this term have 
exhibited solutions based on approximating the 
temperature distribution by a series of Bessel func- 
tions. In the present paper we propose to obtain 
a solution by the method of least squares. The 
boundary value problem will be solved by numerical 
integration of the differential equation employing 
the Runge-Kutta method. This procedure avoids 
(he use of special functions but employs the basic 
concepts necessary to solve the boundary value 
problem. The results obtained indicate that the 
present technique may be employed in other prob- 
lems of a similar nature. 

Let us first put \=2vJK and X£ =x so that (1.2) 
becomes 






d-r 2 ) "™+ 



d*0 
: dr 2 



50 j_d^0 
r dr + X 2 d£ 2 



(1-4) 



while the boundary conditions (1.3) remain un- 
changed except that x is replaced by J. We see from 
(1.4) that when X increases the contribution from the 
term d 2 0/d^ will decrease. By analogy with the 
solution for X= oo it is assumed that 

n — }=iti A n exp (ZnZ)y(r,z n ), (1.5) 

where the constants z n are the eigenvalues and the 
functions y(r,z n ) are the eigenfunctions of the 
boundary value problem 



y 



''+ly'+{§-Zn(l-r>)y y =0, (1.6) 
y(r,\,z) = iorr=l. (1.7) 
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Table 1. The eigenvalues — z n 



n 


>.a = 


X=l 


X=10 


X=100 


X=1000 


\=oo 


1 


2.405 


2. 044369 


6. 744051 


7. 306909 


7. 313523 


7. 313587 


2 


5.520 


5. 187345 


30. 76791 


44. 32333 


44. 60653 


44. 6G947 


3 


8.654 


8. 323452 


59. 50344 


111.9926 


113.9009 


113.9210 


4 


11.792 


11.46137 


89. 47665 


208. 3608 


215. 1669 


215. 2447 


5 


14. 931 


14.60050 


119. 9727 


330. 8685 


348. 3678 


348. 564 


6 


18. 071 


17. 74037 


150. 7501 


476. 6471 


513. 4475 


513. 890 


7 


21. 212 


20. 88073 


181. 6979 


642. 7513 


710. 3849 


711.217 


8 


24. 352 


24. 02119 


212. 7574 


826. 6381 


939. 7287 


940. 551 


9 


27.49 


27. 16301 


243. 9199 


1024. 875 


1196.191 


1201. 87 


10 


30.63 


3G. 30613 


275. 1170 


1224. 774 


1464. 309 


1495. 20 


11 


33.78 


33. 44743 










12 


36.92 


36. 50123 










13 


40.06 


40. 06675 











a For X=0, 2„ are the zeros of the Bessel function Jo(z); the values for X= oo were taken from reference [4]. 



In order to have axial symmetry it is required that 
y(r,\,z) be an even function of r. Further, to deter- 
mine the temperature distribution in (1.5), we 
determine solutions of (1.6) which constitute a 
nonorthogonal set of characteristic functions y(r, 
\,z n ). It is noted that since z occurs quadra tically 
in (1.6) there are two sets of characteristic values. 
It turns out that one set is all positive while the 
other is all negative. However, in evaluating the 
temperature distribution, we wish to have a solution 
which is finite for £ infinite so we discard the 
positive set of characteristic values. In order to 
satisfy the condition 6=d at £=0 we formally 
assume the solution in the form 



71 = 1 



(1.8) 



where /(r) = l in the present problem. When 
X=oo, (1.6) reduces to the Graetz-Nusselt equation 
and the functions -y/r(l—r 2 ), y(r,z n ) constitute an 
orthogonal set and the coefficients are determined 
accordingly. In the present case where X is finite 
the orthogonality property no longer holds. How- 
ever, we may still determine the coefficients in the 
sense of least squares. Consequently, we require 



Cr(l-r 2 )< 



fir)- 



-12A n y(r,z 

n=l 



^y 



dr=minimum. 

(1-9) 



The weight function r(l— r 2 ) has been introduced 
so that the results shall be consistent with the 
Graetz-Nusselt (X=°°) case. It should be noted 
that the solution of (1.6) may be expressed in terms 
of the confluent hypergeometric function. However, 
to obtain the required results it would still be 
necessary to consider the problem of evaluating 



these functions. Instead, we propose to solve the 
differential eq (1.6) directly using numerical inte- 
gration. Differentiating (1.9) with respect to A n 
we get the infinite system of equations 



where 



aA=/3, (1.10) 

«=U= r(l— r 2 ) y(r,z m ) y(r,z n )dr I 

P=^=\_£ r(l-r 2 ) y(r,z m ) f(r)dr~j, 
A=[A n ],m,n= 1,2, .... (1.11) 

2. Determination of z n and y(r,z n ) 

The solution of the infinite system is dependent 
on the determination of the characteristic functions 
y(r,z n ) and the characteristic values z n . The power 
series solution of (1.6) is 



y(r,z n ) = ^2b 2k r 2k , 



k=0 



(2.1) 



where 
b =l,b 2 



-z) b 2k - 2 +zb 2k _±=0. (2.2) 



This representation was used to expedite the 
numerical integration of (1.6). The Runge-Kutta 
method [5] was employed starting at r=0.5 with the 
value of y(0.5,z n ) obtained from (2.1) with trial 
values of z„. The results are given in table 1 (see 
figure 1 for the characteristic functions correspond- 
ing to \=10). 
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Figure 1. 



3. Calculation of the Coefficients Aj 

The integrals in (1.11) necessary for the determina- 
tion of the solution A of (1.10) were evaluated with 
Simpson's rule using the characteristic functions 
y(r, z n ) calculated previously. We give the results 
obtained for 3 cases. Since the matrices are sym- 
metric we omit the elements above the main 
diagonal. 



X=l, 



m = n - 



*a=(10*a mn ) = 



10 7 a=(10 7 <* m J: 



10 7 a=(10 7 <O = 



' L018886 


3951938 

849201 

-176352 

63431 


2472921 

606464 

-133686 

X=10, m=n= 

"239335 

43600 

-13791 

A=1000, m=n- 


1804961 
469689 

= 5 

175401 
37374 

= 8 

133932 

66 110280 

-35 38 

157 158 






L338862 






-242794 






77045 






-32588 
" 947286 


1422473. 




42284 


375196 

15 
-19 

11 

-17 

3 

-14 


S82270 

48645 

-16907 

7130 






-18966 






8546 






-4079 
939335 


138567. 




5 








-12 


234420 

26 170471 
-28 39 

19 -43 
-32 33 
-69 120 




6 

-11 
6 




8 
52 


9379 

-18 


6 
6 



81323. 
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Thus for X=l, if we start with the approximations 
4 !) = — 5 and z { 2 2) = — 6 to z 2 the following successive 
approximations are obtained: 



*2 



1 
2 
3 
4 
5 
6 



-6.0 

-5.0 
-5.219196 

-5. 187761 
-5. 187344 
-5. 187345 



The improved value of z t 
interpolation as follows: 



y(i,4 j) ) 

0. 2362897 

-. 6633393X10" 1 

. 1110549X1CT 1 

. 1455328X10" 3 

-.4580215X10" 6 

-. 2379238X10" 7 

was obtained by linear 



«U + 1)_ 

z 2 



Z 2 



y(l,g,)-y(MJ") ... (1Z )- Q 



Z 2 



y(\A j -")-y(XA S) Y 



The emciency of the method was dependent on the 
choice of the initial values of z. Once some of the 
eigenvalues were known, relatively good starting 
values were easily found by extrapolation. In 
performing the integration (1.6) was written in the 
form 



V 



u =- 



[(z 2 /\ 2 -z)+zr 2 }y 



and the interval Ar= 0.002 was ultimately employed. 

It is noted that as X increases, the matrix a tends 
to a diagonal matrix, i.e., the set of characteristic 
functions y(r, z n ) tends to an orthogonal set with 
respect to the weight factor r(l— r 2 ). Finally, 
examination of the integrals which are the elements 
of a and j8 indicates that since /(r) = l, elements of p 
follow as a simple byproduct of the computation of 
the elements of a. 

Solving the system aA=l3 we get the values of 
A n in table 2. 

With the known values of A n we are able to test 
the emciency of the least square approximation. 
The results are summarized in table 3. 



We conclude from the above that the method 
described yields satisfactory results for practical 
purposes. In conclusion we shall give asymptotic 
representations useful in the regions of small and 
large values of X. These serve to confirm the 
numerical results obtained previously. 



Table 3. ^A n y(r,z n )~f(r) = l 
n = l 



r 


X= 


= 1 


X= 


10 


X = 


100 


X= 1,000 


n=5 


w=8 


n=5 


n=8 


ii=5 


n=8 


n=5 


n=$ 





0.981 


0.918 


1.151 


0.905 


1.188 


0.845 


1.188 


0.864 


0.1 


.937 


1.011 


1.047 


1.009 


1.037 


1.042 


1.033 


1.100 


.2 


.909 


1.004 


0.935 


0.997 


0.921 


0.972 


0.922 


0.961 


.3 


1.004 


0.985 


1.001 


.992 


1.042 


1.033 


1.186 


1.178 


.4 


1.100 


1.023 


1.057 


1.010 


1.040 


0.965 


1. 036 


0.963 


.5 


1.053 


0.971 


0.969 


0.972 


0.920 


1.023 


0.917 


1.036 


.6 


0.965 


1.032 


.942 


1.039 


.997 


0.997 


1.004 


0.974 


.7 


.990 


0.969 


1.076 


0.952 


1.139 


.964 


1.142 


.987 


.8 


.984 


1.024 


1.081 


1.052 


0.992 


1.121 


0.983 


1.126 


.9 


.651 


1.004 


0.682 


0.964 


.539 


0.758 


.528 


0.718 


1.0 



























4. Solution of (1.6) for Small Values of X 

To obtain a solution of (1.6) in terms of Bessel 
functions valid for small values of X we put a\=z 
and obtain 



y" + r-y + [a 2 -a\(l-r 2 )]y=0. 



(4.1) 



Now if X->0, (4.1) becomes y // +r~ 1 y / + a 2 y=0. The 
solution of the limiting problem is J (a s r) where J 
is the Bessel function of the first kind of order zero 
and a s =j s are the zeros of J . 



Table 2. Values of A „ 



n 


X= 


= 1 


X= 


10 


X= 


LOO 


X= 1,000 


n=5 


7i = 8 


71=5 


n=8 


n=5 


n=S 


n=5 


71 = 8 


1 


1. 5779 


1. 6003 


1. 5370 


1. 5423 


1. 478413 


1. 479273 


1.476457 


1. 476472 


2 


-0. 9629 


1. 0527 


-0. 9563 


-0. 9778 


-0. 8126151 


-0. 813778 


-0.806192 


-0. 806222 


3 


.5217 


0. 8172 


.7559 


.8062 


. 599728 


. 606701 


. 588885 


. 588874 


4 


-. 4547 


-. 6658 


-.5894 


-.6865 


-.488371 


-.501968 


-. 475988 


-. 475937 


5 


.2995 


.5489 


.4038 


.5805 


.411076 


. 419506 


. 405103 


. 405024 


6 




-. 4462 




-. 4797 




-.390254 




-. 355635 


7 




.3450 




.3765 




. 348450 




. 318681 


8 




-. 2289 




-. 2565 




-.303294 




-. 287201 
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Thus, if we assume 

a s (\)=a 0s +\ai s +\ 2 a 2s + 
and substitute in (4.1) we find 



(4.2) 



vi(r)- 



1 i 



1 r 2 






«ls 



-j^i ) r ^lO»-| 5 ^WD + g JtW 

7 , 59 



3j! 



15j. 90jf + 90jf 



(4.3) 



Omitting further approximations we show the com- 
parison considering the quantities a ls and a 2 «- 



-a,(\) 



Approx. 



2 (X) 



^4p^roa;. 



1 


2. 04437 


2.04417 


5. 187345 


5. 18755 


2 


3.48773 


3.44802 


9. 75388 


9. 7570 


5 


5. 68968 


6. 0350 


20. 4943 


20. 462 


10 


6. 744051 


15.2 


30. 7679 


32. 52 



5. Solution of (1.6) for Large Values of X 

Let us assume that 

y(r)=y +\- 2 yi +\-%+ . . . 

ZsW = — 20s+2lsX~ 2 +22sX" 4 + • • <• 

We then find that y (r) is the solution of the Graetz- 
Nusselt equation and z 0s are the corresponding eigen- 
values which have been listed in table 1. Substitut- 
ing in (1.6) we obtain a series of characteristic values 
from which we find 



-66.9260 

-2897. 8 



1218. 6 
374 X10 3 



A comparison of the results taking into account z x s 
and z 2s follows: 

X -Si(X) -z 2 {\) 

True value Approx. True value Approx. 



10 


6. 744 


6. 766 


30. 8 


53. 


100 


7. 307 


7. 307 


44. 3 


44. 323 


000 


7. 3135 


7.313 


44. 6065 


44. 606 


00 


7. 313587 




44. 60947 
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ance of Miss Irene A. Stegun in the asymptotic 
developments obtained. 
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